#!/usr/bin/perl
# 2014/05/21 Consider distance effect to calculate P-value

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use File::Basename;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [asocciation file] -o [association file with p-value]\n";
}


my %opt;
getopts("i:o:", \%opt);
my $FILE_aso = $opt{i};
my $FILE_out = $opt{o};
my ($name, $dir, $ext) = &fileparse($FILE_out, '\..*');
my $FILE_tmp = $dir . 'tmp1_' . $name  . '.txt';

#---------------------------------------
# read association file to get total numbers
#---------------------------------------
my %data;
my $INTER_CHR = 0;
my $MEMBER_INTER = 0;
my $fh_aso = IO::File->new($FILE_aso) or die "cannot open $FILE_aso: $!";
while($_ = $fh_aso->getline()){
	s/\r?\n//;
	my ($name1, $loc1, $name2, $loc2, $interaction) = split /\t/;
	my ($chr1, $coord1) = split /:/, $loc1;
	my ($start1, $end1) = split /-/, $coord1;
	my ($chr2, $coord2) = split /:/, $loc2;
	my ($start2, $end2) = split /-/, $coord2;

	if($chr1 ne $chr2){
		$INTER_CHR += $interaction;
		$MEMBER_INTER++;
	}else{
		my $distance = int(abs($start1 - $start2) / 1000) * 1000;
		$data{$distance} += $interaction;
	}

}
$fh_aso->close();

#---------------------------------------
# grouping until group member become more than 100 or >3Mb
#---------------------------------------
my %distance2group;
my %number_in_group;
my %member_in_group;
my $count = 0;
my $group_num = 1;
foreach my $key(sort {$a <=> $b} keys %data){
	$count++;
	if($count > 100 and $key < 3000000){
		$count = 0;
		$group_num++;
	}
	$distance2group{$key} = $group_num;
	$number_in_group{$group_num} += $data{$key};
	$member_in_group{$group_num}++;
}
undef %data;


#---------------------------------------
# read asociation file again and calculate P-value
#---------------------------------------
my $fh_aso2 = IO::File->new($FILE_aso) or die "cannot open $FILE_aso: $!";
my $fh_tmp_out = IO::File->new($FILE_tmp, 'w') or die "cannot write $FILE_tmp: $!";
while($_ = $fh_aso2->getline()){
	s/\r?\n//;
	my ($name1, $loc1, $name2, $loc2, $interaction) = split /\t/;
	my ($chr1, $coord1) = split /:/, $loc1;
	my ($start1, $end1) = split /-/, $coord1;
	my ($chr2, $coord2) = split /:/, $loc2;
	my ($start2, $end2) = split /-/, $coord2;

	my $TOTAL_NUMBER = 0;
	my $TOTAL_COMBINATION = 0;
	if($chr1 ne $chr2){
		$TOTAL_NUMBER = $INTER_CHR;
		$TOTAL_COMBINATION = $MEMBER_INTER;
	}else{
		my $distance = int(abs($start1 - $start2) / 1000) * 1000;
		my $group = $distance2group{$distance};
		$TOTAL_NUMBER = $number_in_group{$group};
		$TOTAL_COMBINATION = $member_in_group{$group};
	}
	my $probability = 1 / $TOTAL_COMBINATION;
	$fh_tmp_out->print("$TOTAL_NUMBER\t$interaction\t$probability\n");
}
$fh_aso->close();
$fh_tmp_out->close();

#---------------------------------------
# calculate P-value using R
#---------------------------------------
my $Rcommand = <<"END";
R --vanilla --slave <<'EOF'
data <- read.table("$FILE_tmp", header=F)

my.pval <- function(M, k, p){
	y <- lgamma(M+1) - lgamma(k+1) - lgamma(M-k+1) + k * log(p) + (M-k) * log(1-p)
	return( exp(y) )
}

line.pval <- function(x){
	pval <- 0
	M <- as.numeric(x[1])
	k <- as.numeric(x[2])
	p <- as.numeric(x[3])
	for(i in k:M){
		pval <- pval + my.pval(M, i, p)
	}
	return(pval)
}

p <- apply(data, 1, line.pval)

asodata <- read.table("$FILE_aso", header=F)
output <- data.frame(asodata, Pval=p)
write.table(output, file="$FILE_out", quote=F, sep="\t", row.names=F, col.names=F)
EOF
END

(system $Rcommand) == 0 or die "R calculation fail : $!";


